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A  sensor  constellation  capable  of  determining  the  location  and  detailed  concentration 
distribution  of  chemical  warfare  agent  simulant  clouds  has  been  developed  and 
demonstrated  on  government  test  ranges.  The  constellation  is  based  on  the  use  of  standoff 
passive  multispectral  infrared  imaging  sensors  to  make  column  density  measurements 
through  the  chemical  cloud  from  two  or  more  locations  around  its  periphery.  A  Computed 
Tomography  (CT)  inversion  method  is  employed  to  produce  a  three  dimensional 
concentration  profile  of  the  cloud  from  the  two-dimensional  line  density  measurements.  In 
this  paper,  we  discuss  the  theoretical  basis  of  the  approach  and  present  results  of  recent 
field  experiments  where  controlled  releases  of  chemical  warfare  agent  simulants  were 
simultaneously  viewed  by  three  chemical  imaging  sensors.  Systematic  investigations  of  the 
algorithm  using  synthetic  data  indicate  that  for  complex  functions,  3-D  reconstruction 
errors  are  less  than  20%  even  in  the  case  of  a  limited  3-sensor  measurement  network.  Field 
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data  results  demonstrate  the  eapability  of  the  eonstellation  to  determine  3-D  eoneentration 
profiles  that  aeeount  for  ~<86%>  of  the  total  known  mass  of  material  released. 


OCIS  codes:  280.0280,  280.4991 

1.  Introduction 

The  improvement  of  ground  test  validation  eapability  in  the  form  of  rapid  and  aceurate 
evaluation  of  sensor  teehnologies  has  reeently  been  targeted  by  the  Joint  Program  Executive 
Office  for  Chemical  and  Biological  Defense  (JPEO-CBD)  as  an  important  technical  objective 
[1].  The  goal  of  the  effort  reported  in  this  paper  is  to  develop  a  system  that  is  capable  of 
providing  accurate  three  dimensional  concentration  distributions  of  chemical  warfare  agent 
simulant  clouds  released  on  test  ranges.  Such  a  system  enhances  the  ability  to  validate  new 
chemical  and  biological  sensing  technology,  as  well  as  provides  a  better  understanding  of  the 
fate  and  transport  of  chemical  releases  through  chemical  reaction  models.  Currently  employed 
technology  serving  this  purpose  makes  use  of  single  pixel  passive  Eourier  Transform  Infrared 
(ETIR)  spectrometers  which  scan  a  viewing  area  every  15  seconds  [2].  Information  on  the 
column  density  of  the  chemical  cloud  obtained  from  the  ETIR  units  is  coupled  to  data  provided 
by  point  sensors  located  across  the  test  area.  The  majority  of  this  information  is  not  currently 
available  for  near-real  time  processing  and  is  collected  manually  after  a  test  is  completed.  One 
goal  of  this  effort  is  the  development  of  capability  to  provide  the  concentration  profiles  in  real 
time  as  well  as  specific  data  presentation  approaches  to  enable  rapid  evaluation  of  the 
performance  of  systems  under  test. 
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The  method  of  eomputed  tomography  (CT)  allows  the  three  dimensional  reconstruction 
of  a  chemical  cloud  density  distribution  from  a  series  of  two  dimensional  measurements  of  the 
line  of  sight  column  density  through  the  cloud  at  different  viewing  angles.  The  use  of  CT  with 
multispectral  passive  infrared  measurements  of  individual  components  in  a  chemical  cloud 
parallels  the  approach  used  in  X-ray  based  medical  imaging  and  can  be  an  important  tool  in 
understanding  chemical  cloud  dynamics. 

In  this  paper,  we  describe  the  development  and  application  of  an  imaging  sensor 
constellation  for  tomographic  chemical  cloud  mapping  that  provides  near-real  time  detection, 
tracking  and  full  3-D  concentration  information  on  chemical  clouds,  both  in  test  situations  as 
well  as  potentially  during  emergency  events.  We  utilize  the  results  of  chemical  imaging 
experiments  conducted  during  the  summer  of  2006  employing  three  LWIR  imaging  Fabry-Perot 
spectrometers  developed  by  Physical  Sciences  Inc.  (PSI)  to  evaluate  the  performance  of  the 
system.  The  instrument  is  known  as  an  Adaptive  InfraRed  /maging  Spectroradiometer,  AIRIS^^ 
(U.S.  Patent  5,461,477).  The  experiments  involved  imaging  a  series  of  controlled  chemical  vapor 
releases  at  Dugway  Proving  Ground  (DPG).  An  illustration  of  the  system  configuration  is  shown 
in  Figure  1.  The  field  study  was  focused  on  the  dissemination  of  several  chemical  warfare  agent 
simulants  viewed  from  three  fixed  locations  approximately  2.8  km  away  from  the  release  point. 

Passive  sensing  of  chemical  vapor  clouds  relies  on  both  the  spectral  signatures  of  the 
target  species  as  well  as  the  radiance  contrast  between  the  vapor  and  the  background  scene.  The 
AIRIS  sensor  units  are  comprised  of  a  long  wave  infrared  (LWIR)  focal  plane  array-based 
camera  which  views  the  far  field  through  a  low-order,  tunable  Fabry-Perot  etalon  [3,4].  The 
tunable  etalon  provides  the  spectral  resolution  necessary  to  resolve  structured  absorption  and 
emission  from  molecular  vapors.  The  focal  plane  array  (FPA)  enables  radiance  measurements  of 
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sufficient  sensitivity  that  the  speeies-speeifie  eolumn  densities  of  ehemieal  elouds  may  be 
determined  with  only  a  eouple  of  degrees  effeetive  temperature  differenee  between  the  vapor  and 
the  baekground.  We  present  a  synopsis  of  the  development  of  the  PSI  tomographie  system  by 
first  deseribing  results  obtained  using  synthetie  data.  This  exereise  has  allowed  us  to  better 
understand  the  absolute  eapability  of  the  algorithm  in  the  ease  of  limited  measurement  network 
(only  three  sensors).  Finally,  we  present  the  results  obtained  from  the  tomographie  analysis  of 
the  infrared  imagery  from  the  DPG  disseminations. 

2.  Experimental  Configuration  and  Theoretical  Basis  of  the 
Approach 

2,1  AIRIS  Sensor  Technology 

The  AIRIS-Wide  Area  Deteetor  (WAD)  instrument  has  been  deseribed  in  detail 
elsewhere  [5,6].  The  basie  AIRIS  optieal  eonfiguration  is  shown  in  Figure  2.  The  eoneept  is 
based  on  the  insertion  of  a  tunable  Fabry-Perot  interferometer  (etalon)  into  the  field-of-view  of 
an  infrared  foeal  plane  array  (FPA).  The  IR  FPA  views  the  far  field  through  the  piezoeleetrie- 
aetuated  etalon  plaeed  in  an  afoeal  region  of  the  optieal  train.  The  tunable  etalon  is  operated  in 
low  order  (mirror  spaeing  eomparable  to  the  wavelength  of  the  light  transmitted)  and  funetions 
as  an  interferenee  filter  whieh  seleets  the  wavelength  viewed  by  the  FPA.  The  optieal 
eonfiguration  depieted  in  Figure  2  affords  a  wide  field  of  view,  high  optieal  throughput,  and 
broad  wavelength  eoverage  at  high  speetral  resolution.  The  interferometer  tuning  time  is  ms 
between  transmission  wavelengths.  Fore-op  ties  and  integrated  blaekbody  ealibration  sourees,  not 
depieted  in  Figure  2,  enable  eontrol  of  the  sensor’s  field-of-regard  and  absolute  radiometrie 
ealibration  of  the  data.  The  AIRIS  units  operated  in  the  study  provided  256  pixel  x  256  pixel 
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imagery  at  2.2  x  2.2  mrad  IFOV  per  pixel.  This  configuration  resulted  in  a  32  degree  x  32  degree 
field  of  regard. 

The  three  AIRIS  units  were  configured  to  acquire  imagery  between  7.9  and  1 1.2  pm  at 
-O.OS  pm  (8  cm"')  spectral  resolution.  The  AIRIS  FPA  acquires  frames  at  200  Hz  and  typically 
8  frames  were  averaged  at  each  observation  wavelength.  The  etalon  tuning  time  to  observation 
wavelengths  is  much  less  than  the  FPA  framing  period  and  the  typical  datacube  acquisition  time 
was  -0.8  seconds.  Datacubes  were  acquired  once  every  -1.7  seconds.  The  noise  equivalent 
spectral  radiance  (NESR)  of  the  data  was  -1.5  pW  cm'  sr'  pm'  (pFlicks)  over  the  entire 
operating  range  and  is  equivalent  to  a  noise-equivalent  temperature  difference  -0.1  K. 
Radiometric  calibration  of  the  data  was  accomplished  using  a  two  point  gain  and  offset 
correction  using  internal  blackbody  calibration  sources  set  for  two  temperatures:  Tcoid  ~  288  K, 
and  Thot  ~  308  K.  A  calibration  of  the  system  was  executed  once  every  30  minutes  for  accurate 
radiometric  measurements.  The  spectral  radiance  of  ambient  temperature  objects  was  typically 
-800  pW  cm'  sr'  pm'  over  the  typical  field  temperature  range  and  the  change  in  radiance  with 
temperature  is  -15  pW  cm'^  sf'  pm''  K''.  Each  256  pixel  x  256  pixel  x  20  wavelength  AIRIS 
datacube  is  used  to  create  a  256  pixel  x  256  pixel  array  of  chemical  vapor  column  density  values 
using  the  approach  described  below. 

2,2  Radiative  Transfer  and  Sensor  Signal  Model 

The  detection  approach  is  based  on  the  change  in  passive  infrared  radiation  received  by 
the  spectrally  resolving  sensor  due  to  the  presence  of  a  chemical  cloud.  The  basic  process  can  be 
described  by  a  three  layer  model,  as  shown  by  Elanigan  [7]  and  is  illustrated  in  figure  3.  The 
total  infrared  radiance  incident  on  the  sensor  at  a  given  wavelength  is  the  sum  of  the 
contributions  from  each  layer  and  is  given  by: 
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where  NbSb  is  the  Planek  radianee  of  the  baekground,  Nc  is  the  Planek  funetion  radianee  at  the 
thermodynamie  temperature  of  the  ehemieal  eloud,  and  Na  is  the  Planek  funetion  radianee  at  the 
thermodynamie  temperature  of  the  atmosphere.  The  quantities  tc  and  tA  are  the  speetral 
transmission  of  the  eloud  and  the  atmosphere  between  the  eloud  and  the  sensor,  respeetively.  The 
first  term  in  Eq.  (1)  deseribes  the  radianee  from  the  baekground  as  attenuated  by  the  ehemieal 
eloud  and  intervening  atmosphere.  The  seeond  term  deseribes  the  radianee  of  the  ehemieal 
speeies  in  the  eloud  as  attenuated  by  the  atmosphere  between  the  eloud  and  the  sensor.  Finally, 
the  third  term  in  Eq.  (1)  deseribes  the  radianee  of  the  atmosphere  between  the  eloud  and  the 
sensor.  The  transmission  of  the  eloud,  layer  2,  is  eomputed  from  the  speetral  properties  of  the 
ehemieal  speeies  eontained  therein; 

tc(^)  =  exp[-^a,(2.)A] 

A  = 

where  Ci  is  the  average  eoneentration  of  the  ehemieal  eompound  over  the  line  of  sight  through 
the  eloud,  I  (/?,  is  the  ehemieal  eolumn  density)  and  ai(k)  is  the  wavelength-dependent 
absorption  eoeffieient.  The  sum  over  index  i  in  Eq.  (2)  is  over  all  speetrally  relevant  ehemieal 
speeies. 

The  differential  radianee  observed  by  the  sensor  as  a  result  of  the  presenee  of  the 
ehemieal  eloud  ean  be  approximated  as: 

m{X)  =  [t^ .  A .  (2,  t;  ) .  (2)  + .  [1  -  A  ] .  (A,  )  +  [1  - 1  j .  n,  (A,  r, )]  -  m)  (3) 
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where  N{A,)  is  the  radianee  measured  at  the  sensor  in  the  absenee  of  the  ehemieal  eloud  (tc^  1) 
and  is  defined  as: 


(4) 


In  order  to  determine  the  ehemieal  simulant  speetral  transmission,  tc  from  measured  sensor 
radiometries,  and  therefore  obtain  its  eolurnn  density,  we  make  use  of  a  fundamental  assumption: 
we  assume  that  the  ehemieal  eloud  attains  atmospherie  temperature  due  to  rapid  mixing  with  the 
surrounding  air.  This  assumption  allows  for  the  approximation  that  Nc~Na  and  Na  can  be 
determined  from  a  measurement  of  the  loeal  air  temperature.  Furthermore,  to  first  order,  this 
assumption  also  relaxes  the  need  to  have  implieit  knowledge  of  the  atmospherie  attenuation,  tA 
between  the  eloud  and  the  sensor.  This  approaeh  has  been  used  in  previous  work  [8].  We  note 
that  atmospheric  temperature  generally  deereases  with  inereasing  height  above  ground  level  and 
that  for  elouds  whieh  rise  high  (several  hundred  meters)  above  the  ground,  it  is  neeessary  to 
address  the  deerease  in  air  temperature  with  altitude  as  the  eloud  mass  will  be  underestimated  if 
the  eloud  temperature  is  presumed  to  be  equal  to  the  ground-level  atmospherie  temperature.  The 
majority  of  the  eoneentration  profiles  estimated  in  this  work  are  ealeulated  for  elevations  within 
tens  of  meters  of  the  ground  however,  so  we  do  not  address  variation  of  air  temperature  with 
altitude  here. 

With  the  assumption  that  the  eloud  temperature  is  equal  to  the  ambient  air  temperature, 
the  eloud  transmission  ean  be  determined  from  Eq.  (1)  and  Eq.  (4)  as  follows: 


"  N{X)-NM) 


(5) 
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where  Nsensor  is  the  sensor  measurement,  and  Na  is  the  Planck  radiance  for  a  blackbody  at 
ambient  air  temperature,  which  can  be  measured  locally.  The  temperature  values  are  typically 
obtained  from  a  weather  station  which  is  positioned  in  proximity  to  the  chemical  release. 

The  most  complex  problem  left  to  be  addressed  in  order  to  solve  Eq.  (5)  is  the 
determination  ofiV(/l) ,  the  scene  radiance  in  the  absence  of  the  chemical  cloud.  As  a  result  of 
the  fact  that  we  do  not  have  a  priori  knowledge  of  the  background  characteristics  (spectral 
features  and  temperature,  Ti),  and  as  a  result  of  the  constant  change  of  the  background 
temperature  with  time,  we  utilize  statistical  methods  for  the  determination  of  iV(/l) .  This 
approach  allows  for  the  identification  of  detector  pixels  containing  the  chemical  simulant,  as  well 
as  the  calculation  of  the  radiance  in  the  absence  of  the  chemical  cloud  even  when  the  simulant 
cloud  is  present  in  the  field  of  view.  The  mathematical  underpinnings  of  this  analysis  method  are 
described  below. 

The  Adaptive  Cosine  Estimator  (ACE)-based  approach  to  target  detection  is  described  in 
more  detail  elsewhere  [9].  ACE  follows  from  an  “unstructured  background”  model.  The 
unstructured  background  model  assumes  that,  while  one  can  determine  an  average  spectrum  for  a 
background  type,  there  is  inherent  spectral  variability  within  the  material  type  and  that,  with  the 
exception  of  lucky  coincidences,  one  cannot  accurately  represent  a  measured  spectrum  by  a 
linear  combination  of  the  average  spectra  of  each  material  type  (as  is  commonly  used  in 
Subspace-Generalized  Eikelihood  Ratio  Test  approach  [9]).  Eor  the  field  data  that  we  have 
investigated  to  date,  the  unstructured  background  model  appears  to  provide  a  better  estimate  of 
the  observed  backgrounds  than  Subspace-GERT  approaches.  Within  this  model,  the  spectra  of 
background  pixels  are  described  as  random  vectors  whose  likelihood  of  occurrence  is  described 


8 


by  a  probability  distribution  function.  The  unstructured  background  model  can  be  described  by 
the  following  expression: 


x  =  ju  +  Sa  +  v  (6) 

where  p,  is  the  mean  spectrum  of  the  spectral/spatial  datacube,  S  is  a  reference  spectrum 
(absorption  coefficients)  for  the  chemical  vapor  of  interest,  and  v  is  a  random  vector  which  is 
characterized  by  a  probability  distribution  function.  The  parameter  ,  as  defined  in  Eq.  (6),  is 
the  differential  radiance  spectrum  due  to  the  presence  of  the  chemical  cloud  as  expressed  in 
Eq.  (3).  In  the  event  that  v  is  normally-distributed  or  elliptically-distributed,  Eq.  (6)  can  be 
expressed  in  terms  of  a  noise-whitened  variable: 

z  =  T;"\x-  ^)  =  T;"^Sa  +  n  (7) 

where  E  is  the  covariance  matrix  of  the  non-target-containing  pixels  and  n  is  a  random  vector 

whose  elements,  when  averaged  over  the  entire  datacube,  have  zero  mean  and  unit  standard 
deviation.  By  working  in  a  noise-whitened  space,  bands  with  low  signal-to-noise  are  de-weighted 
in  the  analysis  and  bands  with  high  signal-to-noise  receive  proportionally  higher  weight, 
hollowing  Eq.  (7),  the  ACE  detection  metric  is: 

T 

Z  Z 

Values  of  Dace  range  from  zero  to  unity,  where  unity  indicates  perfect  correlation  between  the 
test  spectrum  and  the  target  spectrum  and  zero  indicates  no  correlation.  “Target  present”  is 
declared  when  Dace  exceeds  a  user-specified  threshold.  One  of  the  practical  challenges  in 


^ACE  id) 
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applying  ACE  is  determining  whieh  pixels  in  the  seene  to  use  to  ealeulate  p  and  E  .  Eq.  (7)  may 

be  inverted  to  ealeulate  target  speetrum  amplitudes.  Pixels  exhibiting  anomalously  large  values 
of  Mahalanobis  distanee,  typieally  values  greater  than  that  whieh  would  lead  to  inelusion  of  99% 
of  data  from  a  multivariate  normal  distribution,  are  identified  and  the  mean  and  eovarianee  are 
reealeulated  with  those  points  excluded.  In  practice,  when  the  cloud  occupies  a  small  fraction  of 
the  scene,  nominally  a  several  percent  of  the  pixels  or  less,  this  method  of  estimating  the 
background  mean  and  covariance  has  the  effect  of  excluding  most  of  the  target-signature- 
containing  pixels. 

The  least-squares-optimum  estimate  of  target  signature  amplitude  is: 

d  =  (s^r'sy{syr\x-^)  (9) 

Based  on  Eq.  (9),  the  differential  radiance  spectrum  resulting  from  the  presence  of  the 
chemical  in  the  field  of  view  of  the  sensor  ( AN  =  Sd )  is  determined  for  each  datacube 
acquisition  and  for  each  pixel  satisfying  the  Dace  user-specified  threshold.  Consequently,  the 
radiance  at  the  sensor  in  the  absence  of  the  chemical  cloud  is  calculated  as  follows: 

(10) 

where  Nsensor(k)  is  the  measured  radiance  spectrum  and  Sd  is  the  best  fit  to  the  model  as  defined 
above.  Equation  (10)  can  be  used  directly  in  conjunction  with  Eq.  (5)  to  calculate  the  cloud 
transmission  for  each  pixel  identified  as  containing  the  chemical  simulant.  Einally,  following 
Eq.  (2),  the  chemical  simulant  column  density  {p. )  of  chemical  species  i  is  calculated  as: 
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(11) 


,  ^  ,  {dLldT) 

where  (dLjdT')  is  the  derivative  Planek  funetion  with  respeet  to  temperature  averaged  over  the 
wavelengths  in  the  dataeube  and  (N^  -  -  Sd))  is  the  average  differential  radianee 

between  a  blaekbody  at  the  air  temperature  and  the  estimated  baekground  speetrum. 

The  eapability  of  the  deteetion  approaeh  to  aeeurately  measure  the  ehemieal  eolumn 
density  was  initially  validated  under  eontrolled  eonditions  in  the  laboratory.  The  laboratory 
testing  involved  observing  varying  eoneentrations  of  a  ehemieal  in  an  absorption  eell  of  known 
length  viewed  against  a  blaekbody  of  known  temperature.  The  eell  was  filled  with  a  given 
mixture  ratio  of  1,1,1,2-tetrafiuoroethane  (R134a)/nitrogen,  sueh  that  eolumn  densities  of 
100-500  mg/m  were  aehieved.  Directly  behind  the  absorption  cell,  a  large  surface  area 
blaekbody  was  heated  in  order  to  achieve  a  temperature  differential  of 3  °C  with  respect  to  the 
ambient  temperature.  Data  was  acquired  for  several  R134a  column  densities,  and  the  imagery 
was  processed  according  to  the  principles  described  above.  Figure  4  illustrates  R134a 
identification  for  all  pixels  inside  the  absorption  cell.  Figure  5  shows  the  comparison  between  the 
known  column  density  and  the  calculated  chemical  column  density  based  on  the  application  of 
Eq.  (11).  Results  indicate  that  errors  associated  with  ACE-based  column  density  estimation  are 
less  than  10%  when  the  uncertainty  in  the  ambient  temperature  measurement  is  ~  ±  0.3°C. 

The  primary  error  contributions  to  the  chemical  cloud  density  measurement  are  defined 
through  a  simple  error  propagation  analysis  of  Eq.  (11).  The  uncertainty  in  the  column  density 
value  is  evaluated  as  follows: 
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where  Oa  is  the  uneertainty  in  the  ACE-based  statistieal  estimate  of  the  target  speetrum 
amplitude,  ONsensor  is  the  error  assoeiated  with  the  sensor  radiometrie  measurement  (results  from 
the  uneertainty  assoeiated  with  the  radiometrie  ealibration  of  the  sensor  as  well  as  the  sensor 
Noise  Equivalent  Speetral  Radianee  (NESR)),  and  ota  is  the  uneertainty  in  the  ambient 
temperature  measurement.  Eigure  6  illustrates  the  overall  error  eontributions  from  eaeh  of  these 
elements  under  typieal  measurement  eonditions.  In  a  field  deployment  situation,  where  a  single 
weather  station  ean  be  positioned  some  distanee  away  from  the  aetual  chemieal  eloud,  the 
uneertainty  in  the  ambient  temperature  measurement  ean  easily  approaeh  ±  1.0°C.  Consequently, 
it  ean  be  observed  from  Eigure  6  that  the  uneertainty  assoeiated  with  the  temperature 
measurement  introduees  the  largest  error  (~17  pElieks)  to  the  overall  sensor  radiometrie 
measurements. 

3.  Computed  Tomography  Algorithm  for  use  with  Passive  Standoff 
Detection  of  Chemical  Clouds 

Combinations  of  an  integrating  measurement  teehnique  and  CT  have  been  used  in  many 
different  fields  for  the  study  of  diverse  phenomena  [10].  Several  examples  inelude  measurement 
of  soot  volume  fraetions  in  a  flame  based  on  integrated  light  absorption  measurements  [11], 
measurement  of  plasma  emission  intensity  distributions  by  use  of  integrated  emission 
measurements  [12],  and  measurement  of  density  distributions  in  fluid  flows  by  use  of 
interferometry  via  integration  of  the  refractive  index  [13]. 

The  most  difficult  problem  for  optical  tomography  as  applied  to  chemical  cloud 
concentration  mapping  is  the  lack  of  sufficient  projection  data  for  reconstruction;  it  is  not 
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feasible  (cost  prohibitive)  to  acquire  column  density  from  the  cloud  using  very  many  sensors. 
The  ability  to  provide  a  good  tomographic  reconstruction  of  a  chemical  cloud  is  dependent  on 
the  number  of  sensors  used  as  well  as  the  geometric  positioning  of  the  sensors.  The  quality  of  the 
reconstruction  as  a  function  of  these  parameters  was  a  primary  focus  in  determining  the 
tomographic  method  used  for  the  system. 

We  are  given  a  number  of  optical  density  projections  \|/p(x,z)  =  p  (\|/p(x,z)),  where  \|/p  is 
the  column  density  measurement  at  each  pixel  identified  as  containing  the  simulant  in  the  2-D 
plane  orthogonal  to  the  sensor  line  of  sight,  and  from  which  the  chemical  cloud  concentration 
function  C(x,y,z)  is  to  be  reconstructed; 

L 

'^p{x,z)  =  ^C{x,y,z)dy  (13) 

0 

where  the  integration  is  along  the  line  of  sight  y.  Therefore,  given  \|/p  for  p=l,....,P  one  needs  to 
find  an  estimate  of  C(x,y,z).  The  full  3-D  C(x,y,z)  is  generated  by  stacking  each  of  the 
reconstructed  2-D  concentration  distributions,  C(x,y).  Many  algorithms  have  been  developed  for 
the  reconstruction  process,  and  the  choice  of  usage  is  governed  by  factors  such  as  computational 
power  and  speed,  as  well  as  the  spacing  and  projection  data  available.  The  choice  of  CT 
algorithm  for  our  application  was  based  on  the  need  to  generate  the  highest  fidelity 
reconstructions  with  limited  projection  data  and  without  any  a  priori  knowledge  of  the  chemical 
cloud  concentration  distribution.  There  are  generally  two  types  of  CT  algorithms  that  have  been 
extensively  used  in  the  physical  sciences:  transform  based  and  series-expansion  based. 

The  transform  algorithms  are  based  on  the  Fourier  and  Radon  transforms.  The  Fourier 
algorithms  make  use  of  the  projection  (or  central-section)  theorem  [14],  which  relates  the 
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one-dimensional  Fourier  transform  of  a  projeetion  of  an  objeet  to  the  two-dimensional  transform 
of  the  objeet  itself  In  a  limited  projeetion  data  ease  however,  there  are  gaps  in  the  transform 
information  that  must  be  filled  before  the  inverse  transformation  ean  provide  a  reasonable 
reeonstruction  of  the  original  objeet  [10].  The  Radon  transform-based  reeonstruetion  algorithms 
are  iterative  methods  similar  to  the  Fourier  approaches  but  with  the  iteration  carried  out  between 
the  object  domain  and  the  Radon  (projection)  domain.  Support  and  range  constraints  are  applied 
on  the  object  in  the  spatial  domain,  and  the  measured  projections  are  injected  at  each  iteration 
into  the  projection  space  [10].  The  transformations  are  carried  out  by  projection  and  filtered 
backprojection.  The  fundamental  problem  with  transform-based  algorithms  is  that  they  implicitly 
assume  that  the  number  of  projections  is  very  large.  Feng  et  al.  [15]  proposed  an  optimization 
approach  to  try  to  do  a  reconstruction  from  only  3  separate  projection  views.  This  approach 
however,  works  only  for  continuous  and  smooth  object  density  functions.  Furthermore,  the 
algorithms  require  the  line  integrals  to  be  taken  along  uniformly  spaced  and  parallel  projections. 
Even  though  transform-based  algorithms  are  computationally  very  fast,  their  limitations  in 
dealing  with  reduced  projection  data  fields  do  not  make  them  good  candidates  for  use  in  the 
application  proposed  in  this  paper. 

A  better  solution  to  the  reconstruction  of  chemical  clouds  interrogated  by  only  a  few 
sensors  is  accomplished  via  the  use  of  expansion-based  algorithms.  In  the  series-expansion 
formalism,  C(x,y,z)  is  modeled  as  a  linear  combination  of  object  coefficients,  each  of  which 
corresponds  to  the  value  of  the  chemical  cloud  concentration  at  a  point  on  a  rectangular  grid. 
Expansion  algorithms  are  iterative  methods  that  divide  the  volume  of  interest  into  voxels,  and 
one  tries  to  assign  a  concentration  value  to  each  voxel  in  such  a  way  that  the  inferred  path 
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integrals  match  the  observed  integrals.  The  mathematical  definition  of  the  algorithm  is  as  follows 


[10]; 


MNQ 

=  bp(x-Xj,y-yj,z-Zj)ds  (14) 

./  5 

where  Oj  is  the  j*  object  concentration  value  corresponding  to  position  (xj,  yj,  zj)  in  the 
reconstruction  space  and  b  is  a  basis  function  that  weights  the  contribution  of  each  voxel.  The 
(xj_  yj,  Zj)  positions  form  a  rectangular  space  with  M(x  direction),  N  (y  direction),  Q  (z-direction) 
sides.  The  (x,  z)  plane  is  the  plane  of  the  column  density  image  and  the  integral  is  along  y. 
Eq.  (14)  can  be  written  in  simple  matrix  form,  \|/  =  WO; 


Where :  'P  =  Measurement  Vector 
W  =  Projection  Matrix 
O  =  Object  Vector 


(15) 


where  Oj  is  the  local  value  of  the  object  field  C(x,y,z)  at  the  j*  voxel,  and  Wij  is  the  weighting 
factor  which  relates  the  proportion  of  the  j**'  voxel  being  interrogated  by  the  finite  width.  The 
projection  matrix  is  determined  by  the  number  of  sensors  and  the  sensor  positioning  geometry 
used  in  the  measurement.  The  column  density  inferred  from  a  single  sensor  measurement  along  a 
viewing  projection  is  thus  defined  as  a  summation  of  the  concentration  values  for  all  voxels  that 
are  intersected  by  that  projection  weighted  by  the  fraction  of  each  voxel  intersected  by  an 
individual  sensor  pixel  /FOV.  While  the  sensors  provide  fan  beam  projections  through  the 
volume  of  interest,  because  the  cloud  is  in  the  far  field  and  we  are  only  interested  in  estimating 
the  concentration  field  over  a  narrow  range  of  elevation  angles  near  horizontal,  we  approximate 
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the  3-D  distribution  as  a  stack  of  non-interacting  2-D  horizontal  layers  and  the  CT  algorithm  is 
applied  separately  to  each  layer  in  order  to  simplify  the  calculation.  Within  each  2-D  horizontal 
layer  however,  each  viewing  projection  incorporates  a  fan  beam  approach.  This  procedure  is 
illustrated  in  Figure  7. 

In  the  formalism  defined  by  Eq.  (15),  one  would  like  to  solve  the  set  of  P  linear  algebraic 
equations  (one  equation  for  each  measured  line  integral)  in  MNQ  unknowns.  Equation  (15)  may 
suggest  that  in  order  to  determine  concentration  values  (Oj)  one  would  only  need  to  do  a  straight 
forward  matrix  inversion.  However,  one  needs  to  consider  some  important  factors:  1)  the  number 
of  projections  must  be  the  same  as  the  number  of  cells,  2)  for  cases  where  there  are  more 
unknowns  than  equations,  conventional  matrix  inversion  is  no  longer  possible,  and  3)  if  the  data 
suffers  from  the  usual  random  fluctuations  and  noise,  an  exact  solution  will  not  be  possible. 
Clearly,  in  the  case  of  limited  measurement  network  (limited  number  of  sensors),  the  system 
defined  by  Eq.  (15)  is  severely  under-defined.  There  are  several  methods  that  can  be  applied  in 
order  to  provide  the  best  solution  to  such  a  system. 

In  the  Algebraic  Reconstruction  Technique  (ART),  the  solution  to  the  matrix  equation  is 
established  through  an  iterative  process,  in  which,  an  initial  estimate  of  C(x,y)  for  each  cell  in  a 
single  plane  is  given,  followed  by  a  calculation  of  the  column  density,  'Pp,  which  is  compared  to 
the  measured  value  for  the  projection.  A  correction  to  the  initial  value  of  C(x,y)  is  then 
performed,  and  the  process  is  repeated  until  a  given  convergence  criteria  is  achieved.  The  ART 
technique  is  a  widely  used  reconstruction  method  and  is  generally  simple  and  flexible. 
Verhoeven  [10]  presented  a  study  on  the  effectiveness  of  this  algorithm  and  looked  at  its 
reconstruction  ability  by  investigating  the  reconstruction  of  three  simulated  objects  from  three 
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different  geometries.  His  conclusion  was  that  the  best  reconstruction  is  produced  in  the  geometry 
including  fewer  views,  but  which  are  spread  over  a  larger  angle. 

There  are  several  other  expansion-based  algorithms  that  attempt  to  solve  the  matrix 
equation.  According  to  work  by  Todd  et  al.  [16],  reconstructions  performed  by  Maximum 
Likelihood  Expectation  Maximization  (MLEM)  technique  are  found  to  be  statistically  better  than 
ART.  MEEM  is  a  simultaneous  method.  At  each  iteration,  all  voxels  are  corrected 
simultaneously  after  reading  an  entire  set  of  beam  projection  data.  The  tomography 
reconstruction  algorithm  that  we  implemented  for  use  with  IR  passive  standoff  detection  of 
chemical  clouds  is  a  combination  of  the  ART  and  MLEM  techniques.  The  formalism  follows 
Eq.  (15)  and  the  correction  (update)  to  C(x,y)  in  each  iteration  is  determined  by  maximizing  the 
log-likelihood  function  of  observing  a  given  line  integral  through  the  cloud.  The  mathematical 
underpinnings  of  the  approach  are  described  below. 

Eet  us  define  C(xj_  yj,  zj)  =  Oj,  as  the  object  concentration  value  at  the  j*  voxel,  as 
suggested  previously.  In  this  case,  for  each  voxel,  there  is  a  Poisson  distributed  random  variable 
Oj,  with  mean  Oj  that  can  be  generated  independently  (While  the  presumption  that  projection 

measurements  follow  Poisson  statistics  is  not  rigorously  correct  for  the  column  densities 
estimated  here,  in  practice  it  appears  to  be  a  reasonable  assumption  and  it  has  the  benefit  of 
constraining  estimated  concentration  values  to  be  greater  than  zero.): 

PiO.)  =  e-^^^-  (17) 

Ojl 

The  j**'  voxel  produces  a  line  integral  at  \|/p,  which  is  an  independent  Poisson  variable  with 
an  expectation  value  defined  as  (the  projection  of  the  estimated  optical  density  vector): 
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(18) 


MNQ 

=  H^pPj 

J 

Wpj  is  the  probability  (projection)  matrix  identical  to  the  transfer  matrix  from  voxel]  to  the  set  of 
parallel  projections  p  (defined  in  Eq.  (15)).  Since  Oj  are  independent  variables,  a  linear 
combination  of  these  two  variables  follows  Poisson  statistics.  Therefore,  the  likelihood  of  the 
observed  data  is: 


p(Wp)=Y{^ 


(19) 


In  other  words,  this  is  the  probability  under  the  Poisson  model  that  the  given  line  integrals  'Pp  are 
observed  if  the  true  optical  density  is  Oj. 

The  log-likelihood  function  under  MLEM  formalism  is  then  defined  as: 


\n{p{y/^ ))  =  i¥p)+Y.  =  “Z  ) + Z  y^p  ^'^Pp )  ■  Z  ^'^Pp ) 

p=\  p=\  ¥  p-  p=\  p=\  p=\ 
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"Z  Z  ^pPj +11  y^p  Z  ^’pPj  )  ■  Z  ^"^Pp )  (20) 
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MEEM  tries  to  maximize  the  function  in  Eq.  (20).  In  order  to  do  this  maximization,  one 
looks  at  the  first  derivative  of  the  function  with  respect  to  Oj,  and  solves  the  equation  when  the 
derivative  is  zero.  Thus: 


^  P  MNQ  P  U  P 

-^[HP{¥p)]  =  -Z  Z^P.  +Z^i7^ 

‘^7-  p=l  7=1  p=l 
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Dividing  all  sides  by  -  zs  Wpj  ,  then  multiplying  by  Oj  and  combining  with  Eq.  (18)  one  gets: 

p=\  j=\ 
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Eq.  (22)  defines  the  basis  for  the  tomographic  reconstruction,  y/ ^  is  the  projection  of  the 


estimated  concentration  vector  O  at  a  given  iteration  n.  Eet 
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where  T  is  the  multiplicative  coefficient  that  updates  the  value  of  the  j'*'  voxel  at  the  n+1 
iteration.  The  iterations  are  stopped  when  the  following  convergence  criteria  is  satisfied: 


yn 


where  V"  is  the  object  variance  at  the  n*  iteration.  At  each  iteration,  all 


estimated  concentrations  are  updated  simultaneously.  We  initialize  the  concentration  values  by 
assuming  that  they  are  the  same  at  each  voxel.  We  have  also  tried  using  different  initial  guesses 
with  similar  convergence  results.  The  algorithm  typically  converges  in  ten  to  twenty  iterations. 

In  order  to  better  understand  the  capabilities  of  the  tomography  algorithm  in 
reconstructing  chemical  cloud  concentration  distributions  utilizing  a  limited  number  of  views,  we 
considered  sample  objects  of  increasing  spatial  complexity  as  viewed  by  an  increasing  number  of 
sensors,  and  compared  the  reconstructed  functions  with  the  original  analytical  objects.  We  define 
the  metric  for  the  goodness  of  reconstruction  as  the  nearness  [16]  (the  lower  the  nearness  values, 
the  better  the  reconstruction): 
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0  <  Nearness  - 


<1 


(23) 


J 


Oj*  =  j*  voxel  value  of  original  function 

Oj  =  j*  voxel  value  of  reconstructed  function 

Oavg*  =  original  function  average  value 

The  original  analytical  functions  to  be  reconstructed  are  shown  in  Figure  8.  The  first 
function,  Function  A  is  a  simple  2-D  Cosine  Function  (Fig.  8a).  Function  B  (Fig.  8b)  and 
Function  C  (Fig.  8c)  have  the  2-D  Cosine  Function  as  their  basis  but  offer  increased  complexity 
with  the  addition  of  multiple  peaks.  Function  D  (Fig.  8d)  is  a  combination  of  multiple  2-D 
Gaussians  where  the  symmetry  axis  is  not  normal  to  either  of  the  primary  axes  of  the  grid. 

We  first  evaluated  the  algorithm  for  two  orthogonal  sensor  views  (0  and  90  degrees)  with 
128  projections  per  view  for  a  total  of  256  projections  through  the  object.  Reconstruction  of 
Function  A  yields  a  nearness  value  of  0.08.  The  overall  shape  and  maximum  values  are  well 
reproduced  with  errors  less  than  5%,  even  in  the  limited  case  of  2  sensor  views.  The 
reconstruction  result  is  shown  in  Figure  9a  along  with  the  reconstruction  error  distribution. 
Function  B  is  reconstructed  with  a  nearness  value  of  0.17  and  is  shown  in  Figure  9b.  The  overall 
shape  and  maximum  values  are  still  well  reproduced.  There  is  however,  a  clear  decrease  in 
fidelity  at  higher  spatial  frequency,  resulting  in  reconstruction  errors  as  large  as  '-'10%. 

The  reconstruction  of  Function  C  is  shown  in  Figure  9c.  As  a  result  of  the  increased 
spatial  complexity  of  Function  C,  two  sensor  views  cannot  provide  sufficient  projection  data  to 
enable  accurate  estimation  of  the  third  strongest  peak  of  the  function.  It  is  apparent  that  some  of 
the  intensity  is  distributed  along  the  entire  estimation  grid.  The  estimated  distribution  has  a 
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nearness  value  of  0.38  and  estimated  eoneentrations  are  in  error  by  as  muoh  as  40%.  The  need 
for  additional  projection  data  becomes  even  more  evident  when  the  CT  algorithm  is  used  to 
estimate  Function  D.  In  the  case  of  only  two  sensors,  when  cloud’s  axis  of  symmetry  is  not 
parallel  to  one  of  the  viewing  axes,  non-zero  column  density  measurements  can  only  be 
interpreted  as  increased  dispersion.  This  effect  can  be  observed  in  the  reconstruction  of  Function 
D,  as  shown  in  Figure  9d;  the  estimated  distribution  is  characterized  by  an  increased  width  along 
one  of  the  axes.  The  nearness  value  is  0.70  and  is  indicative  of  a  low  fidelity  estimation; 
estimation  errors  are  as  high  as  50%.  These  errors  point  out  the  need  for  an  additional  instrument 
to  reduce  the  ambiguity  in  the  two  sensor  data  set. 

For  comparison  with  the  two  sensor  simulation  we  added  a  third  sensor  whose  line-of- 
sight  was  oriented  at  45  degrees  with  respect  to  the  first  two  sensors  and  used  the  CT  algorithm 
to  estimate  Function  D  again.  The  estimated  function  is  illustrated  in  Figure  10.  The  addition  of 
the  third  sensor  decreases  (improves)  the  nearness  value  from  0.70  to  0.27.  The  errors  associated 
with  the  reconstruction  were  reduced  from  -'50%  to  <20%  and  the  directionality  of  Function  D 
with  respect  to  the  two  orthogonal  primary  grid  axes  is  now  accurately  estimated.  As  might  be 
expected,  the  addition  of  a  fourth  sensor  (at  135  degrees)  reduces  estimation  errors  to  <  10%  and 
improves  the  nearness  value  to  0.10.  The  reconstruction  of  Function  D  with  4  sensors  is  shown  in 
Figure  11. 

In  summary,  our  simulations  lead  us  to  conclude  that  the  tomography  algorithm  cannot 
accurately  reproduce  spatially  complex  clouds  when  only  two  orthogonal  views  are  used  for  the 
cloud  reconstruction.  However,  for  clouds  approximately  described  by  smooth  varying  functions 
with  propagation  along  one  of  the  viewing  centerlines,  estimated  synthetic  cloud  concentrations 
were  accurate  to  within  10%  or  less,  even  when  only  2  projections  views  are  used.  The  addition 
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of  a  third  view  significantly  increases  the  fidelity  of  the  reconstruction  for  spatially  complex 
clouds  with  changing  direction  of  propagation,  resulting  in  errors  of  less  20%.  From  the 
perspective  of  a  perfect  synthetic  projection  data  field,  simulations  indicate  that  4  sensors  may  be 
sufficient  to  provide  full  3-D  reconstructions  of  chemical  cloud  concentrations  that  are  accurate 
to  within  10%,  even  for  very  complex  distributions.  It  is  important  to  note  however,  that  in  real 
field  deployment  scenarios,  the  ability  of  the  sensors  to  achieve  a  highly  accurate  determination 
of  the  cloud  column  density  is  dependent  on  the  level  of  air  temperature  to  background  thermal 
contrast  encountered  by  the  sensors  from  each  viewing  location.  As  a  result,  it  is  very  likely  that 
some  sensors  may  not  provide  enough  column  density  (projection)  data  needed  to  generate 
accurate  3-D  concentration  distributions. 

A  4-sensor  constellation  may  be  degraded  to  a  3-  or  even  a  2-  sensor  constellation 
depending  on  the  quality  of  the  column  density  measurement  associated  with  each  sensor.  As  a 
consequence,  from  a  real  world  application  perspective,  it  is  important  to  provide  a  sufficient 
number  of  sensors  in  order  to  obtain  quality  column  density  data  from  at  least  3-4  viewing 
projection  angles  that  are  needed  for  accurate  reconstructions. 

4.  Dugway  Proving  Ground  (DPG)  Field  Data  and  Tomographic 
Results 

The  data  acquired  by  the  sensor  constellation  during  the  field  trials  at  DPG  in  2006  was 
based  on  observing  explosive  disseminations  of  glacial  acetic  acid  (AA),  tri-ethyl  phosphate 
(TEP)  and  R134a.  All  three  chemicals  display  unique  spectral  signatures  in  the  8-11  pm  spectral 
region  and  are  amenable  to  standoff  detection  by  passive  infrared  sensors  such  as  AIRIS.  The 
majority  of  the  releases  involved  dissemination  of  between  15  kg  and  120  kg  of  material 
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however  three  disseminations  exeeeded  220  kg.  We  present  analysis  of  seven  dissemination 
experiments  here. 

Three  sensors  were  available  for  field  measurements  and  they  were  deployed  in  the  0,  45, 
90  degree  eonfiguration  deseribed  in  the  previous  seetion.  While  in  prineiple  deployment  at  0, 
60,  and  120  degrees  eould  enable  more  aeeurate  estimation  of  the  ehemieal  eoneentration 
distribution,  the  0,  45,  90  degree  eonfiguration  had  the  logistieal  advantage  that  if  one  of  the 
sensors  failed  then  two  sensors  eould  be  immediately  deployed  at  the  0  and  90  degree  loeations 
without  having  to  resurvey  the  deployment  area,  as  would  be  neeessary  if  moving  a  sensor  from 
the  60  deg  or  120  deg  to  90  deg  deployment  position.  The  sensor  eonstellation  was  arranged  sueh 
that  the  field  of  view  of  eaeh  sensor  was  eentered  in  azimuth  on  the  loeation  of  the  explosive 
releases.  Eaeh  AIRIS-WAD  sensor  unit  was  loeated  2.8  km  from  grid  eenter  and  the  geometry 
defined  a  0,  45,  90  degrees  eonfiguration  as  developed  and  validated  using  the  simulation  study 
deseribed  in  the  previous  seetion.  Eaeh  pixel  in  the  image  mapped  to  a  plane  ~6.0  m  x  6.0  m  at 
grid  eenter.  The  absolute  position  of  each  sensor  was  established  to  centimeter-level  accuracy 
using  differential  GPS  measurements.  Extreme  care  was  also  taken  to  align  each  sensor’s  field  of 
view  with  fiduciary  markers  on  the  test  range  to  ensure  that  the  sensor  views  were  co-planar  and 
that  the  projections  through  the  reconstruction  volume  were  known  to  high  accuracy.  Projection 
angles  and  co-planarity  are  believed  to  be  accurate  to  better  than  1  mrad,  i.e.,  significantly  less 
than  the  single  pixel  lEOV,  2.2  mrad. 

Eigure  12  shows  an  example  of  data  products  generated  using  data  acquired 
simultaneously  by  the  three  sensors.  The  cloud  detected  in  Eigure  12a  (top)  originated  from  a 
120  kg  burst  release  of  acetic  acid  and  the  imagery  was  acquired  1  minute  after  the  burst. 
Acetic  acid  detections  correspond  to  the  white  pixels  near  the  center  of  each  image.  The 
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detection  in  Figure  12b  (bottom)  corresponds  to  a  90  kg  TEP  dissemination.  The  Adaptive 
Cosine  Estimator  (ACE),  as  described  in  Section  2,  identifies  and  flags  pixels  which  meet  a 
given  correlation  criteria  between  the  differential  radiance  spectrum  and  the  reference  spectrum 
of  the  species  in  the  cloud.  The  horizontal  lines  in  Eigure  12  determine  the  vertical  overlap 
boundaries  within  which  the  tomography  algorithm  can  be  applied  in  order  to  generate 
concentration  distribution  for  each  altitude  bin  (row  of  pixels)  within  the  detection  grid.  The 
detection  grid  is  defined  as  1.2  km  x  1.2  km  truth  box  in  which  all  three  sensors  have 
overlapping  fields-of-view.  The  full  3-D  cloud  reconstructions  are  obtained  by  stacking  the 
concentration  distributions  from  each  altitude  bin.  Eigure  13  illustrates  the  number  of  detected 
pixels  from  each  sensor  as  a  function  of  time  after  the  initial  release  for  both  the  AA  and  TEP 
disseminations  shown  in  Eigure  12.  The  data  shows  that  the  sensor  constellation  detected  the 
chemical  simulants  within  -'30  seconds  of  the  release.  In  the  case  of  the  AA  release,  the  cloud 
moved  out  of  the  sensors’  fields-of-view  -2.5  minutes  after  the  dissemination  and  re-entered  the 
field-of-view  of  the  sensor  at  the  NE  position  -80  seconds  later.  Typical  sensor  constellation 
detection  times  were  approximately  5-6  minutes  before  the  clouds  cleared  the  detection  grid. 
However,  low  wind  conditions  allowed  some  chemical  releases,  such  as  the  TEP  case  in 
Eigure  12b,  to  be  detected  for  as  long  as  14  minutes.  In  Eigure  12,  the  AA  cloud  was  observed  as 
high  as  100  meters  above  the  ground,  but  the  TEP  cloud  climbed  as  high  as  -300  meters.  (As 
noted  in  Section  2,  concentration  values  estimated  at  elevations  hundreds  of  meters  above  ground 
level  are  likely  underestimated.  Calculations  suggest  that  (effective  thermal  contrast)  may 

be  overestimated  by  as  much  as  10-20%  at  the  highest  cloud  elevations  in  Eigure  12b  if  one 
assumes  that  the  cloud  temperature  is  equal  to  the  air  temperature  at  ground  level;  estimated 
concentrations  would  be  underestimated  by  a  corresponding  amount.)  Most  of  the  detected 
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releases  stayed  relatively  low  to  the  ground  due  to  the  formation  of  low  atmospherie  inversion 
layers  that  are  prevalent  at  DPG  during  the  night.  The  background  scenes  observed  by  the 
sensors  during  the  DPG  releases  typically  provided  air  temperature  to  background  radiance 
thermal  contrast  that  was  >  2  K.  Even  in  cases  when  1  K  <  AT  <  2  K,  the  sensors  demonstrated 
the  capability  to  detect  the  chemical  simulant  clouds. 

Based  on  detection  maps  obtained  from  each  of  the  sensors,  the  chemical  cloud  column 
densities  were  calculated  for  each  pixel  detected  using  the  methodology  previously  described. 
The  column  densities  were  reported  in  units  of  mg/m  .  Column  density  data  was  then  used  as 
input  for  the  tomography  algorithm  in  order  to  generate  concentration  values  in  units  of  mg/m 
associated  with  each  voxel  in  the  3-D  grid.  Furthermore,  the  uncertainties  associated  with  each 
column  density  measurement  were  also  determined  and  propagated  through  the  tomography 
algorithm  in  order  to  generate  the  magnitude  of  the  concentration  error.  As  a  consequence  of  the 
fact  that  the  tomography  algorithm  is  iterative  and  computationally  intensive,  and  since  the 
analysis  described  in  this  paper  was  performed  on  a  slow  IDE  platform,  we  binned  the  sensor 
imagery  from  256  x  256  down  to  128  x  128  in  order  to  reduce  the  computational  time.  As  a 
result,  all  of  the  chemical  cloud  3-D  reconstruction  had  a  nominal  voxel  volume  spatial 
resolution  of  12  x  12  x  12  meters.  Given  the  truth  box  resolution,  and  taking  advantage  of 
observed  low  wind  speeds  (<  2  m/s)  and  slow  moving  clouds,  we  co-averaged  datacubes  over  a 
transit  time  of  ~  20  seconds  prior  to  calculation  of  column  density  maps.  This  approach 
improved  overall  detection  statistics  and  increased  sensor  tomography  overlap.  We  have  recently 
implemented  the  tomography  algorithm  on  an  MKL  C  based  platform  utilizing  a  sparse  matrix 
format,  which  is  currently  capable  of  processing  a  single  256  x  256  slice  in  ~  0.3  seconds,  which 
represents  10  times  improvement  in  computational  speed.  This  recent  accomplishment  allows 
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the  generation  of  tomographie  results  in  near-real  time.  More  reeent  deployment  of  the  real-time 
system  has  demonstrated  the  ability  to  temporally  mateh  dataeubes  within  a  5  seeond  window 
and  generate  3-D  eoneentration  distributions  at  a  rate  of  0.1  Hz.  The  proeessing  period  ineludes 
time  for  all  three  sensors  to  aequire  dataeubes,  reduetion  of  dataeubes  to  d  values  and 
ealeulation  of  p  estimates  at  eaeh  sensor,  wireless  transmission  of  p  values  to  a  eentral 
proeessing  eomputer,  and  the  CT  ealeulation  at  the  eentral  proeessing  eomputer. 

For  demonstration  purposes,  we  focus  only  on  the  acetic  acid  release  shown  in 
Figure  12a,  and  then  provide  a  summary  of  the  complete  analysis.  An  example  of  the  AA 
concentration  distribution  within  the  first  12  meters  above  the  ground  level  is  shown  in 
Figure  14.  Figure  14  is  a  5  x  4  km  representation  of  the  area  indicating  the  position  of  each  of  the 
sensors  along  with  the  1 .2  km  square  detection  grid  in  which  all  three  sensors  have  overlapping 
fields  of  view.  The  calculated  acetic  acid  concentration  distribution  is  overlaid  on  top  of  the  grid 
display.  The  maximum  concentration  detected  is  175  mg/m  .  Figure  15  shows  the 
concentration  distribution  within  the  next  12  m  layer  above  the  ground.  The  maximum 
concentration  is  90  mg/m  ,  indicating  that,  in  the  case  of  the  AA  release,  most  of  the  chemical 
cloud  stayed  relatively  low  to  the  ground.  Figure  16  shows  a  magnified  view  of  the  AA 
concentration  distributions  at  the  two  elevation  levels  above  the  ground. 

In  order  to  validate  the  results  associated  with  the  chemical  cloud  concentration 
distributions,  we  conducted  a  mass  conservation  analysis  for  each  individual  release.  The  only 
known  piece  of  information  regarding  the  dissemination  of  the  simulants  is  the  amount  of 
material  released.  Unfortunately,  the  amount  of  material  present  in  the  gas  phase  was  not  known 
precisely.  While  R134a  is  a  gas  at  ambient  temperature  and  1  atm  pressure,  both  acetic  acid  and 
TEP  have  relatively  low  vapor  pressures  at  ambient  temperature,  '-'16  Torr  and  -'70  mTorr, 
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respectively  and  it  was  not  possible  to  measure  the  partitioning  of  material  between  the  gas  and 
liquid  phase  following  dissemination.  The  amount  of  liquid  used  in  the  dissemination  is  an  upper 
bound  on  the  amount  of  material  which  could  be  present  in  the  gas  phase.  The  tomographic 
analysis  of  chemical  releases  produces  concentration  values  (mg/m  )  for  each  voxel  in  the  3-D 
detection  grid.  Based  on  the  spatial  resolution  of  the  sensor  (dictated  by  the  32  x  32  degrees  FOV 
and  the  2.8  km  standoff  distance),  the  volume  of  each  voxel  was  calculated.  Information  on  the 
voxel  concentration  in  conjunction  with  the  associated  volume  provides  knowledge  of  the  mass 
of  the  simulant  detected  in  each  voxel.  The  total  mass  in  the  3-D  volume  was  then  obtained  by 
summing  over  all  voxels  detected  as  containing  the  simulant.  This  measurement  was  then 
compared  to  the  known  release  mass.  Figure  17  illustrates  the  total  mass  measured  based  on  the 
tomographic  analysis  of  the  release  shown  in  Figure  14a  as  a  function  of  time  after  the  burst.  No 
data  is  shown  after  18:02:30  because  the  cloud  left  the  volume  sampled  by  all  three  sensors.  The 
errors  associated  with  each  mass  measurement  were  calculated  by  propagating  the  column 
density  errors  through  the  tomography  algorithm  and  are  also  indicated  in  Figure  17.  The 
maximum  amount  of  acetic  acid  detected  is  ~  107  ±  30  kg.  When  compared  to  the  actual  amount 
disseminated  during  this  release  (120  kg),  this  value  suggests  that  the  sensor  constellation  in 
conjunction  with  the  tomographic  methodology  is  capable  of  accounting  for  90  ±  25  %  of  the 
total  mass  released  on  the  test  grid.  The  lower  mass  estimates  at  early  times  after  the  burst  are  the 
result  of  the  inability  to  accurately  measure  the  full  column  density  as  a  result  of  the  optical 
thickness  of  the  clouds.  The  simulant  column  densities  in  these  clouds  cannot  be  determined 
using  our  established  radiative  transfer  model.  With  time,  the  chemical  cloud  disperses  and 
reaches  the  optically  thin  limit.  Under  these  conditions,  an  accurate  determination  of  the 
chemical  column  density  is  possible.  Eventually,  the  cloud  disperses  to  the  point  where  column 
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densities  drop  below  the  detection  limit  corresponding  to  the  detection  approach  described  in 
Section  2. 

Similar  analyses  were  conducted  for  several  of  the  other  releases  observed  during  the 
2006  DPG  trials.  A  summary  of  the  mass  estimates  (as  percentage  of  the  total  known  mass 
released)  based  on  the  tomographic  reconstruction  of  the  chemical  clouds  for  each  of  the  runs  is 
shown  in  Figure  18.  It  can  be  observed  that  for  cases  in  which  the  amount  of  chemical  simulant 
released  was  <120  kg,  the  tomographic  constellation  is  capable,  on  average,  of  accounting  for 
~  86%  of  the  total  mass  disseminated.  There  were  three  runs  which  were  characterized  by 
releases  greater  than  240  kg  of  material,  which  resulted  in  average  mass  estimates  only  on  the 
order  of  '-'20%.  Further  examination  of  these  releases  indicate  the  formation  of  optically  thick 
clouds  for  a  significant  fraction  of  the  time  in  which  the  chemical  clouds  are  within  the  detection 
truth  box.  Preliminary  re-analysis  of  the  data  indicates  that  correcting  for  optical  thickness 
effects  increases  the  estimated  mass  by  a  factor  of  2  -  3.  Analysis  conducted  so  far  suggests  that 
the  average  measurement  error  corresponding  to  the  mass  estimate  calculations  are  on  the  order 
of  30%  or  less. 

The  2006  DPG  results  demonstrate  the  ability  of  the  system  to  account  for  large  fraction 
of  the  chemical  simulant  mass  released  on  test  ranges  through  CT  reconstruction  of  cloud 
concentration  distributions.  More  recently  however  (Summer  2008),  the  real-time  tomographic 
sensor  constellation  has  been  deployed  in  conjunction  with  25  ion  mobility  based  point  sensors. 
We  are  in  the  process  of  validating  the  concentration  results  obtained  from  the  tomographic 
reconstruction  of  standoff  passive  AIRIS  data  against  the  raw  concentration  measurements 
collected  with  the  point  sensors.  This  approach  will  provide  a  better  understanding  for  the  ability 
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of  the  proposed  CT  method  to  aeeurately  retrieve  cloud  concentrations  observed  from  true 
standoff  ranges. 

5.  Conclusions 

We  have  defined  the  methodology  and  demonstrated  the  ability  to  use  passive  infrared 
multispectral  imaging  to  track  and  quantify  chemical  clouds  via  computed  tomography  (CT). 
The  CT  algorithm  has  been  demonstrated  to  be  capable  of  3-D  reconstruction  of  chemical  clouds 
using  2-D  column  density  data  with  as  few  as  3  sensors.  Errors  associated  with  the  reconstruction 
methods  are  less  than  20%  even  for  complex  chemical  cloud  concentration  distributions.  The 
statistical  method  for  the  calculation  of  chemical  cloud  column  densities  from  sensor  radiometric 
data  has  been  experimentally  demonstrated  to  produce  chemical  cloud  column  density  values 
that  are  accurate  to  within  ~  10%  or  less.  The  primary  source  of  error  in  the  field  measurements 
was  in  the  determination  of  the  local  air  temperature.  The  total  measurement  error  for  field  data 
is  estimated  to  be  ~  30%.  Unfortunately,  lack  of  ground  truth  data  precludes  assessment  of  the 
absolute  accuracy  of  concentration  distributions  estimated  from  field  test  measurements; 
however,  the  mass  closure  calculations  indicate  that  the  amount  of  material  estimated  using  the 
CT  algorithm  is  in  good  agreement  with  the  known  amount  of  material  released.  The  release  of 
large  amounts  of  simulant  into  small  volumes  during  low  inversion  layers  produced  optically 
thick  clouds.  The  simulant  column  densities  in  these  clouds  could  not  be  determined  using  our 
established  radiative  transfer  model.  This  effect  led  to  an  under-estimation  of  the  mass  detected 
by  the  system. 
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Figure  Captions 

1 .  Imaging  sensor  constellation  for  tomographic  chemical  cloud  mapping. 

2.  Technical  concept  for  Fabry-Perot  imaging  spectrometer.  The  interferometer  is  located  in  an 
afocal  region  of  the  optical  train. 

3.  Schematic  diagram  of  three  layer  radiative  transfer  model. 

4.  Chemical  simulant  identification  in  laboratory  absorption  cell. 

5.  ACE  estimated  column  density  vs.  mass  flow  derived  column  density. 

6.  Error  contribution  to  the  sensor  radiometric  measurement;  1  pElux/pElik  =  1  pW/(cm  sr 
pm). 

7.  Graphical  representation  of  the  reconstruction  technique. 

8.  Test  Eunctions  for  validating  the  reconstruction  algorithm.  A:  2-D  Cosine  Eunction.  B;  2-D 
Cosine  Eunction  as  basis  with  one  additional  peak.  C:  2-D  Cosine  Eunction  as  basis  with  two 
additional  peaks.  D:  combination  of  multiple  2-D  Gaussians  with  random  spatial  distribution 
providing  a  propagation  axis  that  is  not  normal  to  either  of  the  primary  axes  of  the  grid. 

9.  Reconstruction  results  using  2  views  (0,  90  degrees)  and  128  projections/view  a)  Eunction  A, 
b)  Eunction  B,  c)  Eunction  C,  d)  Eunction  D 

10.  Reconstruction  of  Eunction  D  with  3  sensor  views  (0,  45,  90  degrees). 

11.  Reconstruction  of  Eunction  D  with  4  sensor  views  (0,  45,  90,  135  degrees). 

12.  Simultaneous  detection  of  a)  Acetic  Acid  and  b)  TEP  from  all  three  AIRIS-WAD  sensors. 

13.  Number  of  detected  pixels  for  each  of  the  sensors  as  a  function  of  time  after  initial  burst;  a) 
120  kg  Acetic  Acid  release  and  b)  90  kg  TEP  release. 

14.  AA  concentration  distribution  at  0-12  m  elevation  above  ground;  Contour  map  (left)  and 
surface  plot  (right). 
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15.  AA  concentration  distribution  at  12-24  m  elevation  above  ground:  Contour  map  (left)  and 
surfaee  plot  (right). 

16.  AA  eoneentration  distributions  at  0  -  12  m  (left)  and  12  -  24  m  (right)  meters  above  ground. 

17.  Caleulated  total  mass  in  the  3-D  volume  as  a  funetion  of  time  after  the  initial  release  of 
120  kg  of  aeetic  aeid. 

18.  Mass  estimation  based  on  the  3-D  tomographie  reeonstruetion  of  ehemieal  elouds  for  several 
runs  (mass  estimation  expressed  as  pereent  of  the  total  known  mass  released  on  the  test  grid). 
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Error  Temp:  Calibration:  Calibration:  Alpha 
sigma=  Sigma  T  =  Sigma  T=  (ACE) 


+/-1  deg  +/-0.5  deg  +■'-0-3  deg 


Figure  6 


Figure  7 
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Figure  9 
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Figure  10 
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Figure  11 
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Number  of  pixels  detected 


Figure  12 
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Figure  13 
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